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Abstract 

This paper constructs a local generalized finite element basis for elliptic 
problems with heterogeneous and highly varying diffusion coefficient. The 
basis functions are solutions to local problems on vertex patches. The er- 
ror of the corresponding generalized finite element method decays exponen- 
tially with respect to the number of layers of elements in the patches. Hence, 
on a uniform mesh of size H, patches of diameter //log(l jH) are sufficient 
to preserve a linear rate of convergence in H without pre-asymptotic effects. 
The analysis does not rely on regularity of the solution or scale separation 
in the coefficient. This result motivates new and justifies old classes of vari- 
ational multiscale methods. 

Keywords finite element method, a priori error estimate, convergence, multi- 
scale method 
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1 Introduction 

This paper considers the numerical solution of second order elliptic problems 
with strongly heterogeneous and highly varying (non-periodic) diffusion coeffi- 
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cient. The heterogeneities and oscillations of the coefficient may appear on sev- 
eral non-separated scales. It is well known that classical polynomial based fi- 
nite element methods perform arbitrarily badly for such problems, see e.g. IH. 
To overcome this lack of performance, many methods that are based on general 
(non-polynomial) ansatz functions have been developed. Early works [3, IJ, that 
essentially apply to one dimensional problems, have been generalized to the multi- 
dimensional case in several ways during the last fifteen years, see e.g. [fT3llT2i r7]|. 
In these methods the problem is split into coarse and (possibly several) fine scales. 
The fine scale effect on the coarse scale is either computed numerically or modeled 
analytically. The resulting modified coarse problem can then be solved numeri- 
cally and its solution contains crucial information from the fine scales. Although 
many of these approaches show promising results in practice, their convergence 
analysis usually assumes certain periodicity and scale separation. 

For problems with general U° coefficient, the paper [2J gives error bounds for 
a generalized finite element method that involves the solutions of local eigenvalue 
problems. The construction in [161121] depends only on the solution of the original 
problem on certain subdomains. However, the size of these subdomains strongly 
depends on the mesh size. This dependence is suboptimal with respect to the 
theoretical statement given in [illj. that is, for any shape regular mesh of size 
H there exist O i^og{l I H)Y^^^ local (non-polynomial) basis functions per nodal 
point such that the error of the corresponding Galerkin solution uh satisfies the 
estimate \\u - uh\\h^(q,) ^ CfH with a constant C/ that depends on / and the global 
bounds of the diffusion coefficient but not on its variations. The derivation in [fTTll 
is not constructive in the sense that it involves the solution of the (global) original 
problem with specific right hand sides. 

In this paper, we show that such a (quasi-)optimal basis can indeed be con- 
structed by solving only local problems on element patches. We use a modified 
nodal basis similar to the one presented in [ITSl and prove that these basis functions 
decay exponentially away from the node they are associated with. This exponen- 
tial decay justifies an approximation of localized patches. 

The precise setting of the paper is as follows. Let Q c W' be a bounded Lips- 
chitz domain with polygonal boundary and let the diffusion matrix A e L°° {o., Rf^^ 
be uniformly elliptic: 
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Given / e L^(fl), we seek u e V := H^i^) such that 

a(u,v):= [ {AVu,Vv)dx= \ fvdx=:F(v) for all v 6 V. (2) 

The bilinear form a is symmetric, coercive, bounded, and hence, (|2]) has a unique 
solution. 

The main result of this paper (cf . Theorem j8j) shows that the error u - 
of the generalized finite element method, which is based on our new (local) basis 
functions mentioned above, is bounded as follows 

\\A''^V(u-u'^^,)\\<CfH; 

H being the mesh size of the underlying coarse finite element mesh and k « 
log(l///) referring to the number of element layers that form the support of the 
localized basis functions. This result establishes the first error bound for a numer- 
ical upscaling procedure beyond periodicity and scale separation. Moreover, this 
result is stable with respect to perturbations arising from the discretization of the 
local problems. These results give a theoretical foundation for numerous previ- 
ous experiments where the exponential decay of a similar modified basis has been 
noticed, see e.g. [18]. 

The outline of the paper is as follows. In Section [2} we derive a set of local 
basis functions and define the corresponding multiscale finite element method. 
The error analysis is done in Section |3] Section |4] is devoted to the discretization 
of the local problems. Section [5] presents numerical experiments, and Section |6] 
discusses the application of this theory to state-of-the-art multiscale methods. 



2 Local Basis 

In this section, we design a set of local basis functions for the multiscale prob- 
lem under consideration. The construction is based on a regular (in the sense 
of 111) finite element mesh Th of Q. into closed triangles {d = 2) or tetrahedra 



(d = 3). Subsection 2.1 recalls the classical nodal basis with respect to 7~y and 



demonstrates its lack of approximation properties. Subsection 2.3 defines a modi- 
fied (coefficient dependent) nodal basis and analyzes its approximation properties. 
This basis is then localized in Subsection 12.41 



2.1 Classical Nodal Basis 

Let H . Q. ^ R>o denote the 7~//-piecewise constant mesh size function with 
H\t = diam(r) =: Ht for all T 6 Th- The mesh size may vary in space. In 
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practical applications, the mesh Th (resp. its size H) shall be determined by the 
accuracy which is desired or the computational capacity that is available but not 
by the scales of the coefficient. 

The classical (conforming) Pi finite element space is given by 

5// := {v 6 C'^i^) I Vr 6 Th, v\t is a polynomial of total degree < 1}. Qa) 

Let Vh := 5// n y denote the space of finite element functions that match the ho- 
mogeneous Dirichlet boundary condition. Let N denote the set of interior vertices 
of Th- For every vertex x e N, let e S h denote the corresponding nodal basis 
function (tent function), i.e., 

Aj,(x) = 1 and Aj,(y) = for all y ?t x 6 N. 

These nodal basis functions form a basis of Vh- The availability of such a local 
basis is a key property of any finite element method and ensures that the resulting 
linear system is sparse. 

The (unique) Galerkin approximation uh ^Vh satisfies 

a{uH, v) = F{v) for all v 6 Vh- (|3}b) 

The above method ([3]) is optimal with respect to the energy norm |||-||| := : = 
IIA^^^V • \\ii(Q,) on V which is induced by a, 

\\\u - uh\\\ = min - v//||| . (4) 

Assuming that the solution u is smooth, the combination of (|4]) and standard in- 
terpolation error estimates yields the standard a priori error estimate 

\\\u - Uh\\\ < C||//||L»(n)l|V^M||t2(n). 

This estimate states linear convergence of the classical finite element method ([3]) 
as the maximal mesh width tends to zero. However, the regularity assumption 
is not realistic for the problem class under consideration. Moreover, even if the 
coefficient is smooth, it may oscillate rapidly, say at frequency s'^ for some small 
parameter s- In this case, the asymptotic result is useless because V^m may os- 
cillate at the same scale, a fact that is hidden in the constant ||V^M||i2(j^) ^ - 
Unless H < s, the above finite element space is unable to capture the behavior 
of the solution neither on the microscopic nor on the macroscopic level. In what 
follows, we present a modification of this method that resolves this issue. 
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2.2 Quasi Interpolation 

The key tool in our construction will be some bounded linear surjective (quasi-) 
interpolation operator '■ V ^ Vh- The choice of this operator is not unique and 
a different choice leads to a different multiscale method. We have in mind the fol- 
lowing modification of Clement's interpolation [jlOl which is presented and ana- 
lyzed in [8, Section 6]. Given v eV, 3//v := T,.xeN(-^Hv)(x)Ax defines a (weighted) 
Clement interpolant with nodal values 

(3ffV)(x) := (f^vA.dx) Kf^A.dx) (5) 

for X e N. The nodal values are weighted averages of the function over nodal 
patches tOx '■= supp A^. Since the summation is taken only with respect to interior 
vertices A'^, this operator matches homogeneous Dirichlet boundary conditions. 

Recall the (local) approximation and stability properties of the interpolation 
operator 3// [8, Section 6]: There exists a generic constant C^^ such that for all 
V 6 V and for all T e 7"// it holds 

Ht'Wv - ^HV\\L2iT) + I|V(V - ^HV)\y(T) < C3 J| Vv||i2(,,) (jSJa) 

with cot := U{^ e Th I T n K d}. The constant Ca^ depends on the shape 
regularity parameter p of the finite element mesh Th (see ([13]) below) but not on 
diamT. 

The above interpolation operator is not a projection, i.e., Vh g Vh does not 
equal its interpolation ^hVh in general. However, the particular choice gives rise 
to the following lemma. 

Lemma 1 There exists a generic constant C'~^ which only depends on p but not 
on the local mesh size H, such that for allvu 6 Vh there exists v 6 V such that 

3ff(v) = V// and |||v||| < C-^ |||v//||| and supp v c supp v//. (|6jb) 

Proof. For every nodal basis function Ax, x e N, there is a /3x 6 H^^itOx) such that 
"^Hifix) = Ax and |||;0.v||| < C"^ IPxIII with some constant C'~^ that does not depend 
on X and H. Note that Px may be chosen from the finite element space that is 
related to some uniform refinement of Th- 

Given vh = Z.veyv^ffW^A 6 Vh, v := Vh + Z.veyv(^//W - {^hVh){x))Px 6 V 
has the desired properties. The interpolation and support properties are obvious. 
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The stability follows from 
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supp /l.vl^'' the inverse inequality ||v//-3//V//||^„(7^) < 



and (|6la|). 



where we have used | 

In what follows the interpolation operator ([5]) may be replaced by any lin- 
ear bounded surjective operator that satisfies ( |6[a[ )-( [6[b| ). Hereby, ( |6[b[ ) may be 
relaxed in the sense that suppv is not necessarily a subset of suppv// but that 
supp V \ suppv// covers at most a fixed (small) number of element layers about 
supp Vh. 



2.3 Multiscale Splitting and Modified Nodal Basis 

Let 3// : V ^ V// be a quasi interpolation operator according to the previous 
subsection. Then the kernel of 3// 

y^= {v e y I = 0} 

represents the microscopic features of V, i.e., all features that are not captured by 
Vh- Given v eVu, define gv 6 by 

a(5v, w) = a(v, w) for all w e V^. 

The finescale projection operator g : V// ^ leads to an orthogonal splitting 
with respect to the scalar product a 

V=V^'® where V^' := (Vh - Wh). 

Hence, any function u e V can be decomposed into 6 V^^ and e V\ 
u = u^^ + u\ with a(M™', u^) = 0. Since dim VjJ" = dim Vh, the space V^^ can 
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be regarded as a modified coarse space. The superscript "ms" abbreviates "mul- 
tiscale" and indicates that V^^, in addition, contains fine scale information. The 
corresponding Galerkin approximation e V^^ satisfies 

aiu"^', v) = Fiv) for all v e (7) 

The error (u - m™) of the above method Q is analyzed in Section 13.1 



Finally, we shall introduce a basis of V^'^. The image of the nodal basis func- 
tion Ajc under the fine scale projection 5 is denoted by = 5^.v 6 V\ i.e., (px 
satisfies the corrector problem 

a{(px, w) = a{Ax, w) for all w e V^. (8) 

A basis of V^^ is then given by the modified nodal basis 

{Ax - (Px \xeN}. (9) 

In general, the corrections (f>x of nodal basis functions Ax, x e N, have global 
support, a fact which limits the practical use of the modified basis (|9]) and the 
corresponding method (|7]). 

2.4 Localization 



In Section 3.2 we will show that the correction (px decays exponentially fast away 
from X. Hence, simple truncation of the corrector problems to local patches of 
coarse elements yields localized basis functions with good approximation proper- 
ties. 

Let 6 N. Define nodal patches of k-th order o*^ about x e Nhy 

a)x,i := snppAx = int(U{r eTnlxeT}), 

0Jx,k := int(U {T eTnlTn Ux,k-i ^ 0}) , ^ = 2, 3,4 . . . . 

Define localized finescale spaces V^(a)x,k) := {v e | v\n,\aj,,t = 0}, x e N, 
by intersecting with those functions that vanish outside the patch (jL)x,k- The 
solutions (f)x,k 6 V^(o^x,k) of 

a((f>x,k, w) = a(Ax, w) for all w e V\a)x,k), (1 1) 

are approximations of 0, from ([8]) with local support. 
We define localized multiscale finite element spaces 



V^l = spm{Ax - (Px,k \xeN}^V. ([TIja) 
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The corresponding multiscale approximation of ([2]) reads: find e such 
that 



a(M^^^, v) = F(v) for all v e V^^. ([TIjb) 

Note that dim V™^ = \N\ = dim Vh, i.e., the number of degrees of freedom of 
the proposed method < [T2| ) is the same as for the classical method (|3]). The basis 
functions of the multiscale method have local support. The overlap is proportional 



to the parameter k. The error analysis of Section [3^ suggests to choose k « log ^ 

Remark 2 The localized modified basis fiinctions could be localized fiirther to 
vertex patches co^, x e N, by simply multiplying them with the classical nodal 
basis functions; for any x e N and any y e N C0x,h define (f>l : = Ay(px,h The 
generalized finite element space which is spanned by those O {{[og{\ I H)Y^ local 
basis functions per vertex has similar approximation properties as ^ (see [^5]). 



3 Error Analysis 

This section analyzes the proposed multiscale method in two steps. First, Subsec- 



tion 3.1 presents an error bound for the idealized method (|7]). Then, Subsection 
3.2 bounds the error of truncation to local patches and proves the main result, that 
is, an error bound for the multiscale method ( [T2] ). 

As usual, the error analysis depends on the constant p > which represents 
shape regularity of the finite element mesh 7~y; 

p := maxpr with pj := for T e Th, (13) 

teTh diam T 

where Bj denotes the largest ball contained in T. 
3.1 Discretization Error 

Lemma 3 Let u eV solve ^ and u^'' e V™' solve (|7]). Then it holds 

\\\u-u^%<C'fC:s,a-'\\Hf\\L2^a) 
with constants Coi and Ca^ that only depend on p. 



Proof. Recall the (local) approximation and stability properties ( |6[a[ ) of the inter 
polation operator 3//. Due to the splitting from Section 2.3 



it holds u - = u^. 
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Since 3//m^ = 0, the application of ( |6[a| ) and Young's inequality yield 
|||i<l||' = F(i<0< ll/lb(r)ll"'-3//"'lb(r) 

for any e > 0. Note that there exists a constant Coi > that only depends on p 
such that the number of elements covered by coj is uniformly bounded (w.r.t. T) 
by Col- The choice e = concludes the proof. ■ 

Remark 4 Substituting by the modified Clement interpolation operator pre- 
sented in / fTOl/ allows one to replace the term WHfWiHci) in Lemma ([3]) by data 

oscillations [Yj^eN \\H{f - fxWii^^^ J with some weighted averages fx of f on 
cOx, X 6 N; we refer to [10 Section 2] for details. 



3.2 Error of Localized Multiscale FEM 

First, we estimate the error due to truncation to local patches. We wiU frequently 
make use of cut-off functions on element patches. 

Definition 5 For x e N and d < D eN, let r]i'^ : Q. ^ [0, 1] be a continuous and 
weakly differentiable function such that 

iri'^'')L., = 0, (Uga) 
(rii''')\iiW_..o = h and ([T4}b) 

vr e Th, \\Vn.-c\\L-iT) < Cco(£> - dy'H-' ^c) 

with some constant Qo that only depends on p. For example, one may choose 
Tji''^ e S H with nodal values 

rf^^{x) = forallxeNn to^, 

rii'^ix) =1 forallxeN<^{n\ oj^^d) , and (15) 
yfx'^ix) = m{D - dy^ for all x e N O doJx,d+m, m = 0,1,2, . . . ,D - d. 

Lemma 6 For all x e N, k,£ > 2 eM the estimate 

holds with constants Ci, C2 that only depend on p but not on x, k, €, or H. 
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Proof. Letx e N and{,k>2 e N. Observe that 



-0./^||P<lll<^.v-v||P = |||f,-v|||2^.^^ + |||0 



(16) 



holds for all v 6 V^(a)x,(k) using Galerkin orthogonality. 

Let ■■= 1 - 77^-1^+''^*=-^ with a cutoff function as in DefinitionH 

According to ([6|b]), there exists e V such that ^nibx) = 3//(^.v0.v), lll^vlll ^ 
Il|3//(^x0.r)lll, and suppC/?,,) c Hence, v := ^;,^.v - b^c e V^(co^,ek) and 

^ c^„C3„ (lll0.IIL,,«\<.,,*-iHi + lll^.v'^-IL.,«\-v/,.-i)) • 

Since 3//0x = 0, the upper bound of the interpolation error ( |6[a[ ) and ( |14[c[ ) yield 

<c'^\\W"' 



with C2 := 1 + CoiQqCj^Ca and Ca '■= maxTeTn ^(Xt)- This leads to 

^<^2lll0.IIL.,,rW.,..(.-.K• 
The combination of ([16]), with v = ^^(('x - bx, and ( [T7] ) yields 



0. 



(17) 



(18) 



Further estimation of the right hand side in ( 18 ) is possible using cutoff func- 



(19) 



tions rjj := rj^^ ^ I'+i-^^' (cf. Definition [5J), j = 2, 3, . . . , A; - 1. Observe that 

= J (AV0„V(77ti0.)}dx - 2j^77^.i0,<AV0,,V77^.i)dx 

Let, according to Lemma [IJ b^xk-i) be chosen such that ^nbxxk-i) = ^ninl^i^Px)- 
Then ril iCp^ - b^xk-i) £ ■ Since | supp(V/i^) n supp(?7<:_i)| = and supp(V?7i.„i) = 
ii>x,(k-v)i \ <^x,(k-2)i+i, the first term on the right-hand side of ( [191 ) can be rewritten 
as 

JjAV(Px,V{7]l,(f>x))dx 

= J (AV(l)x,V(r]l_i(t>x-bxxk-i)))dx + J {AV(px,Vbxxk-i))dx 



/ 

Jn 



(20) 



{AV(f)x,'^bxxk-i))dx 



'<^x,(k-i)t\'JJx,{k-2)l'+\ 
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,(t-l)f\cj.v,(*:-2)f+l 



With rfj ■.= \T\ ^ ^ 77^_j Ax we have 



<c 



< C3„ (ll(/7li - vDh-iT) IUMt + l|V(/7Li)||L"(r) ^m\L2(T)) 

< 2C3j|V(77,_i)|L»(r) (diam(r) IUMt + ^\\<Px - ^h{<Px)\\lHt) ■ 



Thus the property ( |14[c[ ) of the cutoff function and the upper bound of the inter- 
polation error (|6[a|) yield 



<^x,(k-l)t\<JJy:.{k-2)t+l 



,(t-2)f)' 



(21) 



where C\ only depends on Cj^, Ceo, Coi, and Ca- The same arguments allow one 



to bound the second term on the right-hand side in ([19]) 
Jn 



< 2 IIV/7^-illL~(r)IIA'^'V0,|L2(r)||A>/V.llL2(r) (22) 



TeTn'- rc(jj,-_(j-_i)f\iL)v(j-_2)f+i 

<c';e-'\\A'i^vuli 



L2(n\(j.,,(*-2y)' 



where C" only depends on C^^, Ceo, and Ca- The combination of ([T9])-(|22]) yields 

ll^^^'V0.,||^,^,^^^^^_^^^^ < Cir^||A>/^V0.||^,^^^^^^^^_^^^^, (23) 

where C\ := C[+ C". For 7 = A; - 2, . . . , 2, a similar argument (with rik-i replaced 
by rjj) yields 



l|l4l/2^ 



-i)f)' 



(24) 



Starting from ( [23] ), the successive application of ( |24] ) for j = - 2, - 3, . . . , 2 
proves 



< (Cir>)*^-2||A>/2y^^ 



(25) 



Combining (]T8]) and (25 1, we finally obtain the assertion. 



Lemma 7 T/zere /5 a constant C3 ?/?a? depends only on p and PI a, but not on \N\, 
k, or t such that 



1 

v(x)(0, - 0,,«) < c^itkf Yj ^'(^) 1 1 k-v - 



x&N 



lk\ 



xeN 
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Proof. For x e N, let = 1 - rf!^^^-^^*^ (cf. Definition H). By Lemma [T] there 
exists a function hy^ e V such that for any w e it holds ^h^x = - 
^x)w), suppibx) c supp(3ff((l - ^x)w)) c a)xjk+3 \ oj^jk, and lll^xllL,,„+3\a;,.„ ^ 
|||3//((1 - &)w)\\L^,,^,\oj^,,. We note that w-^^w-Z?;, 6 with support outside 
a>x,(k , i.e., w-^.^w-Z7;,) = a(Ax, w-^^w-bx) = and vv-^^^w-Z^.v) = 0. 

With w = X.veA^ v(a:)(0.:, - cl)jc,ek) 6 we have 

IlkllP = ^ vix) a{(px - (px,ek^ ^xW + bx) 

xeM 

< J] \vix)\ \\\(f>x - (l>x,ik\\\ ■ (lUxwlw + c:^ |||3,,((i - ^.)w)|IL^„.^3) 

xeN 

xeN 

< 2C'^C^,,liYj\v{x)\ - cPx,a\\\ ■ ll(V^,)(l - 3//)w)||^2(, 

xeN 

+ 2C^^C3„ ^ IvWI 

xeN 

( 

Y^v\x)\\\(Px-(px,. 



:,ft • lllH^II 



a 



■jklW ■ \\M 



xeN 



dl2 



\xeN 



1/2 



\W\\\ , 



where Co^{£kY is an upper bound on the number of patches tOx/k that overlap a 
single element in the mesh. The result follows by dividing by |||w||| on both sides. 



Theorem S Let u e V solve ^ and, given £,k > 2 e N, let m™'^^ e V™'^^ solve 
([12]). Then 

\\\u - < Ca'WHj'U^^a) (Ci/£n\\f\\L2(n) + C:s,a-'\\Hf\y^^y 

holds with C\ from Lemma^and a constant C4 that depends on j3/a, and p but 
not on H, k, €, f, or u. 

Proof. Let m^^^^ := YixeN "™(-^) (^a: - 4>x,ck)^ where xeN, are the coeffi- 

cients in the basis representation of Due to Galerkin orthogonality, Lemma 
|3} Lemma |7j and the triangle inequality, 



< Cj,ff--||H/llJj,n| + C3(«)" "hM^ Ilk' - A. 



2NU .IIP (26) 



xeN 
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The application of Lemma [6] yields 

xeN xeM 

Furthermore, we have 

2 

LHT) 



xeN 



TeT 



xeTnN 



TeT 



TnN 



= ci 



xeN 



2 

mm 



2,.msi|2 I /^r 



xeN 



2 

L2(n) 



where C[^^ and Ca^ depend on p and Cp is the constant from Friedrichs' inequality. 
This yields 



H,tk\ 



\(k-2)l2 |||,,ms 



<CmT\\L-ia)a\lkr\C,l€) 



dl2. 



U 

\(k-2)/2 



H 



(27) 



where C4 only depends on C2, C3, C[^^, Cp, and Ca^. The assertion follows readily 
by combining ( [261 ) ■ 



Remark 9 T/ze error estimate in Theorem^contains a factor ||//~'||L"(n)- How- 
ever, its influence on the total error can be controlled by choosing the localization 
parameter k proportional to log(l/|i//~' For non-uniform meshes, it is rec- 

ommended to vary the choice of the localization parameter in space according to 
k « log ^. 



4 Discretization of the Fine Scale Computations 

In this section focus on how to compute numerical approximations to the local 
basis functions - (px/k and there by to the multiscale solution w^**^^. In order 
to do this, we need to extend the error analysis of Section 3 to a fully discrete 
setting. There is a lot of freedom in choosing different finite elements and different 
refinement strategies, see e.g. [fT5l[T6l . We will focus on a very simple and natural 
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approach, also considered in [[T7l. We assume that the local basis functions are 
computed using subgrids of a fine scale reference mesh, which is a (possibly space 
adaptive) refinement of the coarse grid 7~y. 

More precisely, let Th be the result of one uniform refinement and several 
conforming but possibly non-uniform refinements of the coarse mesh Th- We 
introduce h : Q. ^ R>o as the 71,-piecewise constant mesh width function with 
hi := h\, = diam(0 for all t eTh. We construct the finite element space 

5/, := {v 6 C°(Q) I Vt 6 ThiQ),v\t is a polynomial of total degree < 1}. 
We let 6 V/, := 5/, n Hj^i^) be the reference solution that satisfies 

aiuh, v) = F(v) for all v 6 V,,. (28) 
Locally on each patch we let 

Vl(oj,^k) := V\co,^k) n Vh = {v 6 ni 3^v = and v\n\^^, = 0}. (29) 
The numerical approximation 0^ ^ 6 V;^(<y.v,/t) of the corrector 0^ ^ is determined by 

a(0j.^, w) = a(A^, w) for all w 6 Vl(a)j,^k)- 
We denote the discrete multiscale finite element space 

V^f = span{i, - \xeN}. 
The corresponding discrete multiscale approximation u^^j^ e ^H^t fulfills 

a(u^f, V) = F(v) for all v 6 V'^f. (30) 
Theorem 10 LetueV solve ^ and let u'^f^ e solve Then 

III" - <«||| ^ Cmr'UL-im {Cd£)-^\\f\\L2^a)+C,s,a-\\Hf\y^a)M\\u - uu\\\ . 

Remark 11 (Multiscale splitting by nodal interpolation) Having discretized the 
fine scale computation, i.e., having replaced the infinite dimensional space V by 
some finite element space Vf, c C°(Q) we are allowed to replace the Clement- 
type interpolation by classical nodal interpolation. This leads to the variational 
multiscale method in [16]. Because nodal interpolation satisfies the conditions 
( |6[a[ )-( [6jbl ), Theorem 10 establishes an a priori error bound for the multiscale 
method 4761/. 
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Remark 12 (Estimates for the fine scale error) The finite element space Vh may 
be replaced by any finite element space that contains V/„ e.g., by piecewise poly- 
nomials of higher order. The third part in the error bound in Theorem 10 can 
be bounded in terms of data, mesh parameter h, and polynomial degree using 
standard a priori error estimates. For example, if A 6 W^'°°{Q.) (bounded with 
bounded weak derivative) and s is the smallest present scale, i.e., ||VA||z,»(n) < e~^ 
the third term in the error bound in Theorem^^may be replaced by the worst case 
bound Chs'^ for a first-order ansatz space Vh (see /|20|/]. It is shown in / |2Q|/ that 
for highly varying but smooth coefficient A, higher order ansatz spaces are supe- 
rior 

Remark 13 (Periodic coefficient) Let Q. be some square or cube, f e L}{0), A 
be smooth and periodic, A{x) = A(x/s), with some small scale parameter e > 0, 
and let u^ denote the corresponding solution of ([2]). Choose uniform meshes Th 
and Th with H > s > h and k x \og(H~^). With regard to the previous comment. 
Theorem \W\yields the error bound 

III nis,/i 1 1 1 ^ / T T , /i \ 

With h ~ sH the error in the approximation becomes independent of the fine scale 
oscillations without any so-called resonance effects as they are observed, e.g., in 

m. 

Remark 14 (Solution of the local problems) The local problems need to be solved 
in the spaces This is a standard finite element space with the additional 

constraint that the trial and test functions should have no component in Vu- In 
practice this constraint is realized using Lagrange multipliers. 

The resulting coarse scale system of equations is of the same size as the orig- 
inal problem, dim{V^^^) = dimVn and it is still sparse. The number of non-zero 
entries will be larger and depend on k. Note however that the non-zero entries in 
the stiffness matrix decay exponentially away from the diagonal. 



Proof of Theorem 
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\nh - M^'^'^'lll + ~ "//«||| ^'^d follow the arguments from the proof of Theo- 



We use the triangle inequality |||m - "^^111 - HI" ~ "''HI ^ 



rem [8] simply replacing V by V/, and using Lemmas 15 17 and 18 below (discrete 



versions of Lemmas [3|[6[ and|7]) to bound the last two terms. 

In the remaining part of this Section, A < B abbreviates an inequality A < C B 
with some generic constant < C < oo that does not depend on the mesh sizes H, 
h, and the localization parameters. 
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Lemma 15 (Discrete version of Lemma |3) Let Uh e Vh solve ( [28] ) and i/^''^ e 
y^^'^ solve ([30]) with k large enough so that oj^j, = Qfor all x e N. Then 

|||M/,-<'''|||<c;;'c3„a~'ii///ib(n) 

holds with constants Coi and C^^ that only depend on p. 

Proof. Note that wjj := Uh- u^''^ is the unique element of := n Vi, such that 
a{i/j^,v) = F(v) for all v/ 6 Vl- The Lemma follows from the same arguments in 
the proof of Lemma [3j ■ 

To establish discrete versions of Lemmas [6] and [7] we are facing the technical 
difficulty that the product of v 6 Vh and some cut-off function 77 from Definition |5] 
is not necessarily an element of Vf,. However, the subsequent lemma shows that 
the product rjv can be approximated sufficiently well by elements from Vh. 

Lemma 16 For all x e N, D > d e N, and corresponding cut-off function rji'^ 



defined in ([15]) there exists some v e VI{cOx,d+i) such that 



Furthermore the statement also holds ifrji'^ is replaced by 1- rji'^ and v e V^iQ. \ 

Proof. Let xeN, D>deMhe fixed and define 77 := 7]i'^. Let ^h ■■ V n 
C(Q.) V^ be the nodal interpolant with respect to the mesh 71,. Recall its (local) 
approximation and stability properties 

||V(v - ^hV)\\LHt) ^ hW^MtHt) and W^nvWiHt) ^ MkHt) 

for all polynomials v. According to Lemma [T| there exists some bx 6 V'' such that 
^H{bx) = ^H(Mr](f>% \\\bx\\\ < |||3//(3ft(#J-))|||, and supp(^;,) c aJx,D+i \ (Ox,d-i- 
Hence, v := 3/,(#^) - bx e Vl{a)x,D+i)- Since ^hMx = for = \t\~^ lr](y)dy, 
we get 



2 III .;, ^ , .h^ , ll|2 

2 



||k'-v||| =\\\ncf>':-Mricl>l) + b 



< 

fe71,:fca).v,o\iiJ.,,rf 

teTh:tcwx^D\<^y.d TeTH.TciI>_y^D+\W.xj-\ teTi,:tcT 

teTh-.tccbx^oW^j TeTH.T'Z(b,,D+i\'^,.d-i teTi,:tcT 



^ Yj l|/^ll^(r)IIV/7llz 



iL~(r)|||Vx|||r 
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The first part of the Lemma now follows using the property ( |14[c[ ) of rj. The second 
part concerning l-rj follows using the same argument but with v e V^i^ \ <^x,d-\)- 



Lemma 17 (Discrete version of Lemma |6|) For all x e N, k,£ > 2 e N the 

estimate 



A.' 



Ck\ 



<C2 



k-2 

-11 \\\ M 



e 



holds with constants Ci, C2 that only depend on p but not on x, k, i, h, or H. 



Proof. Let := 1 - with as in equation ^ in Defini- 

tion [s] Then there exists a v 6 Vlico^/i,) such that, 



v\\\ < 

< 



Furthermore, using the same argument as in Lemma |6| 



■^'''<^x/k-l\'^i.((k-l}+i 



■^'''iJJx.nWx.Clk-D+l 



which yields 



Ck\ 



< 



< 



C(k-\)+l 



Now let rij := rf^^ i)+i ''^ (^f. Definition [s]), j = 2,3, . . . ,k - \ and note that 
The second term can be treated exactly as in Lemma [6] and, hence, bounded by 



2 

r2||U':||L .We make use of Lemma 
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to bound the first term. There exists 



'x,i{k-l)+y 



) such that 



a{(pWk-i^\) <\\yi\\\ ^ llkll'^'-HIU^■'IIKII^o^ 

vr.V IK ITA/ \\\^ ^\\\clC{k-\)\'^e{k-2)+\ III III lll^-'^lllf^\tJ.v/(*:-2) 

The final assertion follows by similar arguments as in the proof of Lemma [6j ■ 

Lemma 18 (Discrete version of Lemma |7]) There is a constant C3 depending only 
on p and fi/a, but not on \N\, k, or t such that, 



xeN 



h |||2 

'xlk\\\ 



xeN 
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Proof. For x e N,let = ^ - 
(pxjk)- We have, 



ek+\,ek+2 



(cf. Definition [5| and z = Y^xeN ^(-^) (spx - 
v{x) {(Px - = Yj ^(^) «(<^x - + ( 1 - ^x)z) = 1 + 11. 



xeN 



xeN 



The first term I := YjxeN ^(^) '^(4>x ~ 0v [v ^^z) can be treated in exactly the same 
way as in the proof of Lemma|7} We focus on the second term. Due to Lemma 16 
there exists awe Vl{Q. \ oja) such that 

n := J] v{x) a{<p'l - (1 - - w) 



xeN 



< 



< 



xeM 



2 \\\^h _ .h 

X Yxlk\ 



\xeN 



> 


1/2 / 


||2 




) 


V 


^ 


1/2 / 


||2 




) 


V 



J] 111(1 will 

xeN 

Z 



1/2 



''ijJx.fki-lWfk-H 



J 



\xeN 

The result follows immediately. 



2 III J./! Jr 



x/k\ 



1/2 



Z ■ 



5 Numerical Experiments 

Numerical experiments shall validate our theoretical results from the previous sec- 
tions. 



5.1 Experimental setup 

Let Q be the unit square and the outer force / = 1 in Q. Consider three different 
choices for the scalar coefficient Ai, A2, A3 with increasing difficulty as depicted 
in Figure [T] The coefficient Ai = 1 is constant. The coefficient A2 is piecewise 
constant with respect to a uniform Cartesian grid of width 2'^. The values in 
each grid cell are chosen in the range [1/20,2]; the contrast /3(A2)/a{A2) < 40 
is moderate. The coefficient A3 is piecewise constant with respect to the same 
uniform Cartesian grid of width 2"^. Its values are taken from the data of the 
SPEIO benchmark, see http : //www. spe . org/web/csp/. The contrast for A3 
is large, ;S(A3)/a(A3) « 4 ■ 10^. Consider uniform coarse meshes of size H = 
2 \ 2"^, . . . , 2~^ of Q. as depicted in Figure [2| Note that none of these meshes 
resolves the rough coefficients A2 and A3 appropriately. 
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Figure 1: Scalar coefficient used in the numerical experiment: Ai (left), A2 (mid- 
dle), A3 (right). 



Figure 2: Uniform triangulations of the unit square. 

The reference mesh Th has width h = 2~^. Since no analytical solutions are 
available, the standard finite element approximation ui, 6 Vh on the reference mesh 
Th serves as the reference solution. The approximations are compared with this 
reference solution only. Doing this, we assume that w/, is sufficiently accurate. All 
fine scale computations are performed on subsets of Th. 

5.2 Results for the energy error 

Figure |3] depicts the energy errors of the new multiscale method and the classical 
PIFEM (see (|3])) with respect to the same coarse mesh. Depending on the coarse 
discretization scale H, the localization parameter k is chosen to be r21og(l///)]. 
The logarithmic dependence on l/H is motivated by our a priori analysis. The 
choice of the constant 2 is based on numerical tests. It turns out that, in all experi- 
ments, this choice leads to the desired linear textbook convergence (rate -1/2) of 
the energy error (w.r.t. to the number of degrees of freedom A^dof = W\ ~ H~^) 
related to the sequence of multiscale approximations. Pre- asymptotic effects are 
not observed. In particular, the performance of our method does not seem to be 
affected by the high contrast present in A3. Whether our estimates on the decay of 
the corrector functions are sub-optimal or have worst-case character with respect 
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to contrast is an issue of present research. 

Observe that the classical PIFEM suffers from the lacks of approximability 
and regularity and converges only poorly for the rough coefficients A2 and A3. 

5.3 Results for the error 

Figure [4] shows errors of the new multiscale method and the classical PIFEM. 
Again, the choice of the localization parameter k = r21og(l///)] yields the opti- 
mal convergence rate -1 for our method in all experiments (w.r.t. to the number 
of degrees of freedom Ndof = W\ « H~^) without any pre- asymptotic behavior. 

More importantly, we observe that the error between Uh and con- 
verges nicely at a rate close to -3/4. This is remarkable because, is a 
true coarse approximation. "^hU^^ element of the coarse PI finite element 
space. Hence, it cannot capture microscopic features of the solution. The rate of 
convergence is limited by for some 5 e [0, 1] which is related to the regularity 
of the solution {u e H^^" for some s e [0, 1]). However, in contrast to the classical 
PIFEM approximation uh, '^rU^'^ converges at some rate better than -1/2 and, 
hence, allows to approximate the macroscopic behavior of the solution very accu- 
rately with only very few degrees of freedom. Note that the storage complexity 
of the modified basis is of order OQr^ log 1 ///) whereas its interpolation can be 
stored in 0{H^ log 1 ///). Once the system matrix of the multiscale method is as- 
sembled, '^nu^i' can be computed without using the modified basis whereas this 
would be required to represent u^l\ 

6 Application to Multiscale Methods 

In this section we discuss three multiscale methods and how the presented analysis 
relates to each of them. 



6.1 The Variational Multiscale Method 

The variational multiscale method was first introduced in [fT3]| . The function space 
V is here split into a coarse part (standard finite element space on a coarse mesh), 
in our case Vh, and a fine part, in our case V'. The weak form is then also decou- 
pled into a coarse and a fine part. The method reads: find u eVn and u' 6 such 
that, 

a(u, v) + a{u , v) = F(v) for all v 6 Vh, 
a{u', v') = F{v') - a(u, v') for all v' e V^. 
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Figure 3: Relative energy errors - / \\\uh\\\ (□ solid red) with localiza- 

tion parameter = r21og(l///)] and \\\uh - uh\\\ I \\\uh\\\ (O dotted blue) vs. number 
of degrees of freedom A^dof ~ for different coefficients: A\ (top), A2 (middle), 
A3 (bottom). The dashed black line is N'.^F . 
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Figure 4: Relative errors \\uh - u^''^\\l\\uh\\ (□ solid red), \\uh - ^HU^f\\l\\uh\\ (* 
solid red) with localization parameter = r21og(l///)] and (O dot- 

ted blue) vs. number of degrees of freedom A^dof * for different coefficients: 
Ay (top), A2 (middle), A3 (bottom). The dashed black line is A^^J^. 
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The fine scale solution is further decoupled over the coarse elements T e Th and 
approximated using analytical techniques. Note that the fine scale solution u' is 
an afline map of the coarse scale solution u. If we let u' « Mu + m and plug this 
in to the method we get a coarse stiffness matrix of the form a{v + Mv, vv), i.e., a 
non-symmetric bilinear form for a symmetric problem. 

6.2 The Multiscale Finite Element Method 

In [fT2l the multiscale finite element method was first introduced. Here modified 
multiscale basis functions are computed numerically on sub-grids on each coarse 
element individually. The basis functions fulfill: find (p^j 6 H^{T) 

a(Ax - (p„j, v) = for all v e HliT) and for all T e Th- 

Here homogeneous Dirichlet boundary conditions are used on the boundary of 
each element T, i.e., the local problems are totally decoupled. To get a more 
accurate method one can improve the boundary conditions using information from 
the data A. A larger domain can also be considered (this procedure is referred to 
as over-sampling), see e.g. [fT2l| . Note that since the coarse scale basis functions 
are modified (both trial and test space) the resulting method is symmetric. 

6.3 The Proposed Method 

The modified basis function construction given by equation ([8]j9j) was first in- 
troduced in a variational multiscale framework in [[T4l [TSll . In these papers the 
Scott-Zhang interpolation was used in the analysis and nodal interpolation in the 
discrete setting for the numerical examples. The method has also been used to 
modify basis functions in the spirit of the multiscale finite element method, see 
e.g. ffTBl. The exponential decay in terms of layers of coarse elements have been 
demonstrated numerically in several works, see e.g. [fT6l [T8l [TTll . In all experi- 
ments considered in these papers, the fixed choice of two or three layers of coarse 
elements in the patches appeared sufficiently accurate for engineering purposes. 
The first paper where the exponential decay was analyzed is [17J. The theory of 
iterative solvers was used to show convergence to a discrete reference solution. 

The proposed method has been extended to convection dominated problems 
and problems in mixed form [iTSil as well as time dependent problems [fT9ll . A 
posteriori error bounds have been derived and adaptive algorithms designed where 
the local mesh and patch size are chosen automatically in order to reduce the error. 
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6.4 Application of the Presented Analysis 

The convergence proof in this paper generalizes the results of [fTTll and gives a 
valid bound also as h ^ independent of the patch size and coarse mesh size. 
The proof does not rely on regularity of the solution and gives a very explicit ex- 
pression for the rate of convergence. The present analysis confirms the numerical 
results in [|T6irT8llT7]l and gives together with [17J the method the solid theoretical 
foundation that has previously been missing. The analysis also justifies the use of 
a posteriori error bounds for adaptivity [ [TSlfTSlI because we can now prove that the 
quantities measured on the patch boundary decays exponentially in the number of 
coarse layers. 

For the variational multiscale method this result says that it is important to 
allow larger subgrid patches than just one coarse element. This will result in over- 
lap but the local problems are totally decoupled and we have in previous works 
demonstrated how adaptivity can be used to only solve local problems where it 
is needed, see for instance ifTSlfTSll . For the multiscale finite element method the 
analysis is not directly applicable since the fine scale space is not used. It is the 
decay in this space which we have proven to be exponential (in number of coarse 
layers of elements in the subgrid). If this decay is not present, inhomogeneous 
boundary conditions are instead needed for the subgrid problems. To the best of 
our knowledge, such constructions have only been proved to be accurate in special 
settings, e.g., periodic coeflicients. 
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